Task4

df <- read.csv("winequality-red.csv", header = TRUE, as.is = FALSE)
head(df)

##   fixed.acidity volatile.acidity citric.acid residual.sugar chlorides
## 1           7.4             0.70        0.00            1.9     0.076
## 2           7.8             0.88        0.00            2.6     0.098
## 3           7.8             0.76        0.04            2.3     0.092
## 4          11.2             0.28        0.56            1.9     0.075
## 5           7.4             0.70        0.00            1.9     0.076
## 6           7.4             0.66        0.00            1.8     0.075
##   free.sulfur.dioxide total.sulfur.dioxide density   pH sulphates alcohol
## 1                  11                   34  0.9978 3.51      0.56     9.4
## 2                  25                   67  0.9968 3.20      0.68     9.8
## 3                  15                   54  0.9970 3.26      0.65     9.8
## 4                  17                   60  0.9980 3.16      0.58     9.8
## 5                  11                   34  0.9978 3.51      0.56     9.4
## 6                  13                   40  0.9978 3.51      0.56     9.4
##   quality
## 1       5
## 2       5
## 3       5
## 4       6
## 5       5
## 6       5

Кoрoтко о данных

Это данные о физико-химических свойствах вина, например, о кислотности, содержании диоксида серы, сульфатов и тд, а также качество вина (измеряемое баллами от 0 до 10).

1 - фиксированная кислотность

2 - летучая кислотность

3 - лимонная кислота

4 - остаточный сахар

5 - хлориды

6 - свободный диоксид серы

7 - общий диоксид серы

8 - плотность

9 - рН

10 - сульфаты

11 - алкоголь

Первичный анализ данных

Описательная статистика

summary(df)

##  fixed.acidity   volatile.acidity  citric.acid    residual.sugar  
##  Min.   : 4.60   Min.   :0.1200   Min.   :0.000   Min.   : 0.900  
##  1st Qu.: 7.10   1st Qu.:0.3900   1st Qu.:0.090   1st Qu.: 1.900  
##  Median : 7.90   Median :0.5200   Median :0.260   Median : 2.200  
##  Mean   : 8.32   Mean   :0.5278   Mean   :0.271   Mean   : 2.539  
##  3rd Qu.: 9.20   3rd Qu.:0.6400   3rd Qu.:0.420   3rd Qu.: 2.600  
##  Max.   :15.90   Max.   :1.5800   Max.   :1.000   Max.   :15.500  
##    chlorides       free.sulfur.dioxide total.sulfur.dioxide    density      
##  Min.   :0.01200   Min.   : 1.00       Min.   :  6.00       Min.   :0.9901  
##  1st Qu.:0.07000   1st Qu.: 7.00       1st Qu.: 22.00       1st Qu.:0.9956  
##  Median :0.07900   Median :14.00       Median : 38.00       Median :0.9968  
##  Mean   :0.08747   Mean   :15.87       Mean   : 46.47       Mean   :0.9967  
##  3rd Qu.:0.09000   3rd Qu.:21.00       3rd Qu.: 62.00       3rd Qu.:0.9978  
##  Max.   :0.61100   Max.   :72.00       Max.   :289.00       Max.   :1.0037  
##        pH          sulphates         alcohol         quality     
##  Min.   :2.740   Min.   :0.3300   Min.   : 8.40   Min.   :3.000  
##  1st Qu.:3.210   1st Qu.:0.5500   1st Qu.: 9.50   1st Qu.:5.000  
##  Median :3.310   Median :0.6200   Median :10.20   Median :6.000  
##  Mean   :3.311   Mean   :0.6581   Mean   :10.42   Mean   :5.636  
##  3rd Qu.:3.400   3rd Qu.:0.7300   3rd Qu.:11.10   3rd Qu.:6.000  
##  Max.   :4.010   Max.   :2.0000   Max.   :14.90   Max.   :8.000

Удалим из рассмотрения free.sulfur.dioxide, так как есть total.sulfur.dioxide и удалим citric.acid, так как есть fixed.acidity.

df <- df[-3]
df <- df[-5]

Виды признаков

Таблица внизу показывает, какую моду имеют наши количественные признаки. Отсюда можно сделать вывод о том, являются они непрерывными или нет. Построим matrix plot для того, чтобы увидеть особенности в наших данных.

library(lattice)
library(ggplot2)
library('GGally')
ggpairs(df, title="correlogram", columns=c(1:9), upper = list(continuous = "points"), diag = list(continuous = "barDiag"))

Прологарифмируем 1, 3, 4, 5, 8, 9.

dfl <- transform(df, fixed.acidity=log(fixed.acidity), residual.sugar=log(residual.sugar), chlorides=log(chlorides), total.sulfur.dioxide=log(total.sulfur.dioxide), sulphates=log(sulphates), alcohol=log(alcohol))
ggpairs(dfl, title="correlogram", columns=c(1:9), upper = list(continuous = "points"), diag = list(continuous = "barDiag"))

Удалим единичные outliers.

dflo <- dfl
dflo[rownames(dflo)[dflo$volatile.acidity > 1.5 | dflo$chlorides < -4 | dflo$total.sulfur.dioxide > 5.1 | dflo$alcohol > 2.7],] <- NA
ggpairs(na.omit(dflo), title="correlogram", columns=c(1:9), upper = list(continuous = "points"), diag = list(continuous = "barDiag"))

names(dflo)[names(dflo) == 'fixed.acidity'] <- 'log_fixed.acidity'
names(dflo)[names(dflo) == 'residual.sugar'] <- 'log_residual.sugar'
names(dflo)[names(dflo) == 'chlorides'] <- 'log_chlorides'
names(dflo)[names(dflo) == 'total.sulfur.dioxide'] <- 'log_total.sulfur.dioxide'
names(dflo)[names(dflo) == 'sulphates'] <- 'log_sulphates'
names(dflo)[names(dflo) == 'alcohol'] <- 'log_alcohol'

Нормальность данных:

library(ggpubr)
ggqqplot(dflo$log_fixed.acidity, ylab = "log_fixed.acidity")

ggqqplot(dflo$volatile.acidity, ylab = "volatile.acidity")

ggqqplot(dflo$log_residual.sugar, ylab = "log_residual.sugar")

ggqqplot(dflo$log_chlorides, ylab = "log_chlorides")

ggqqplot(dflo$log_total.sulfur.dioxide, ylab = "log_total.sulfur.dioxide")

ggqqplot(dflo$density, ylab = "density")

ggqqplot(dflo$pH, ylab = "pH")

ggqqplot(dflo$log_sulphates, ylab = "log_sulphates")

ggqqplot(dflo$log_alcohol, ylab = "log_alcohol")

shapiro.test(dflo$log_fixed.acidity)

## 
##  Shapiro-Wilk normality test
## 
## data:  dflo$log_fixed.acidity
## W = 0.98536, p-value = 1.234e-11

shapiro.test(dflo$volatile.acidity)

## 
##  Shapiro-Wilk normality test
## 
## data:  dflo$volatile.acidity
## W = 0.97928, p-value = 2.158e-14

shapiro.test(dflo$log_residual.sugar)

## 
##  Shapiro-Wilk normality test
## 
## data:  dflo$log_residual.sugar
## W = 0.85681, p-value < 2.2e-16

shapiro.test(dflo$log_chlorides)

## 
##  Shapiro-Wilk normality test
## 
## data:  dflo$log_chlorides
## W = 0.82868, p-value < 2.2e-16

shapiro.test(dflo$log_total.sulfur.dioxide)

## 
##  Shapiro-Wilk normality test
## 
## data:  dflo$log_total.sulfur.dioxide
## W = 0.98794, p-value = 3.061e-10

shapiro.test(dflo$density)

## 
##  Shapiro-Wilk normality test
## 
## data:  dflo$density
## W = 0.99131, p-value = 4.189e-08

shapiro.test(dflo$pH)

## 
##  Shapiro-Wilk normality test
## 
## data:  dflo$pH
## W = 0.99327, p-value = 1.206e-06

shapiro.test(dflo$log_sulphates)

## 
##  Shapiro-Wilk normality test
## 
## data:  dflo$log_sulphates
## W = 0.95749, p-value < 2.2e-16

shapiro.test(dflo$log_alcohol)

## 
##  Shapiro-Wilk normality test
## 
## data:  dflo$log_alcohol
## W = 0.946, p-value < 2.2e-16

Нормальность отсутсвует :( Но объем выборки позволяет говорить об асимптотической сходимости методов.

Посмотрим на корреляцию между признаками.

library(ggcorrplot)
library(ppcor)
dflo <- na.omit(dflo)
corr <- round(cor(dflo[c(1:9)]), 1)
ggcorrplot(corr, method = "circle", lab = TRUE)

ggcorrplot(pcor(dflo[c(1:9)])$estimate, method = "circle", lab = TRUE)

Сильная зависимость наблюдается между log_fixed.acidity и density, а также log_fixed.acidity и pH.

library(lm.beta)
lm.df.oecd <- lm.beta(lm(data = dflo, quality ~ log_fixed.acidity + volatile.acidity + log_residual.sugar + log_chlorides + log_total.sulfur.dioxide + density + pH + log_sulphates + log_alcohol))
summary(lm.df.oecd)

## 
## Call:
## lm(formula = quality ~ log_fixed.acidity + volatile.acidity + 
##     log_residual.sugar + log_chlorides + log_total.sulfur.dioxide + 
##     density + pH + log_sulphates + log_alcohol, data = dflo)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.67653 -0.35971 -0.05412  0.44655  1.91522 
## 
## Coefficients:
##                           Estimate Standardized Std. Error t value Pr(>|t|)    
## (Intercept)               47.95847      0.00000   22.85206   2.099 0.036006 *  
## log_fixed.acidity          0.49971      0.12368    0.21900   2.282 0.022635 *  
## volatile.acidity          -0.89892     -0.19787    0.10357  -8.680  < 2e-16 ***
## log_residual.sugar         0.10168      0.04479    0.06365   1.597 0.110370    
## log_chlorides             -0.22137     -0.08816    0.05872  -3.770 0.000169 ***
## log_total.sulfur.dioxide  -0.07963     -0.06939    0.02490  -3.199 0.001408 ** 
## density                  -48.90313     -0.11373   23.17792  -2.110 0.035024 *  
## pH                        -0.09064     -0.01738    0.19333  -0.469 0.639249    
## log_sulphates              0.80965      0.22449    0.08466   9.563  < 2e-16 ***
## log_alcohol                2.63132      0.32207    0.29004   9.072  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6427 on 1582 degrees of freedom
## Multiple R-squared:  0.3639, Adjusted R-squared:  0.3602 
## F-statistic: 100.5 on 9 and 1582 DF,  p-value: < 2.2e-16

Значимыми являются с ***. согласно модели quality будет больше при меньшем значении 4 признаков (volatile.acidity, log_chlorides, log_total.sulfur.dioxide, density; максимальный по модулю– volatile.acidity) и большем значении 3 (log_sulphates, log_alcoho, log_fixed.acidityl; максимальный по модулю– log_alcohol). Остальные коэффицинты

AIC: Выбор модели для линейной регрессии

Также выбор модели можно осуществлять на основании информационных критериев, которые базируются на функции максимального правдоподобия для оценок параметров и учитывают число параметров.

В случае функции stepAIC, AIC вычисляется по формуле \[\text{AIC} = 2k - 2 \log \mathcal L\], где \(k\) — число параметров модели, \(\mathcal L\) — функция максимального правдоподобия для модели.

“Backward” алгоритм stepAIC будет убирать переменные по одной, пока не будет достигнут минимум функции AIC.

lm.df.oecd.full <- lm.beta(lm(data = dflo, quality ~ log_fixed.acidity + volatile.acidity + log_residual.sugar + log_chlorides + log_total.sulfur.dioxide + density + pH + log_sulphates + log_alcohol))
stepAIC(lm.df.oecd.full, direction = "backward")

## Start:  AIC=-1397.63
## quality ~ log_fixed.acidity + volatile.acidity + log_residual.sugar + 
##     log_chlorides + log_total.sulfur.dioxide + density + pH + 
##     log_sulphates + log_alcohol
## 
##                            Df Sum of Sq    RSS     AIC
## - pH                        1     0.091 653.55 -1399.4
## <none>                                  653.46 -1397.6
## - log_residual.sugar        1     1.054 654.51 -1397.1
## - density                   1     1.839 655.30 -1395.2
## - log_fixed.acidity         1     2.151 655.61 -1394.4
## - log_total.sulfur.dioxide  1     4.226 657.68 -1389.4
## - log_chlorides             1     5.871 659.33 -1385.4
## - volatile.acidity          1    31.119 684.58 -1325.6
## - log_alcohol               1    33.998 687.45 -1318.9
## - log_sulphates             1    37.778 691.23 -1310.2
## 
## Step:  AIC=-1399.41
## quality ~ log_fixed.acidity + volatile.acidity + log_residual.sugar + 
##     log_chlorides + log_total.sulfur.dioxide + density + log_sulphates + 
##     log_alcohol
## 
##                            Df Sum of Sq    RSS     AIC
## <none>                                  653.55 -1399.4
## - log_residual.sugar        1     1.449 655.00 -1397.9
## - density                   1     3.576 657.12 -1392.7
## - log_total.sulfur.dioxide  1     4.144 657.69 -1391.3
## - log_chlorides             1     5.871 659.42 -1387.2
## - log_fixed.acidity         1     7.587 661.13 -1383.0
## - volatile.acidity          1    31.276 684.82 -1327.0
## - log_sulphates             1    38.914 692.46 -1309.3
## - log_alcohol               1    43.791 697.34 -1298.2

## 
## Call:
## lm(formula = quality ~ log_fixed.acidity + volatile.acidity + 
##     log_residual.sugar + log_chlorides + log_total.sulfur.dioxide + 
##     density + log_sulphates + log_alcohol, data = dflo)
## 
## Coefficients:
##              (Intercept)         log_fixed.acidity          volatile.acidity  
##                 54.00449                   0.58040                  -0.90063  
##       log_residual.sugar             log_chlorides  log_total.sulfur.dioxide  
##                  0.11194                  -0.21461                  -0.07845  
##                  density             log_sulphates               log_alcohol  
##                -55.26861                   0.81480                   2.56140

Информационный критерий AIC убирает всего 1 предикатор– pH.

Удалим признак log_fixed.acidity (большая корреляция со многими признаками) и построим модель линейной регрессии ещё раз.

corr <- round(cor(dflo[c(2:9)]), 1)
ggcorrplot(corr, method = "circle", lab = TRUE)

ggcorrplot(pcor(dflo[c(2:9)])$estimate, method = "circle", lab = TRUE)

Теперь не наблюдается сильных корреляций (как обычных, так и частных) между признаками.

lm.df.oecd.full <- lm.beta(lm(data = dflo[c(2:10)], quality ~ volatile.acidity + log_residual.sugar + log_chlorides + log_total.sulfur.dioxide + density + pH + log_sulphates + log_alcohol))
stepAIC(lm.df.oecd.full, direction = "backward")

## Start:  AIC=-1394.4
## quality ~ volatile.acidity + log_residual.sugar + log_chlorides + 
##     log_total.sulfur.dioxide + density + pH + log_sulphates + 
##     log_alcohol
## 
##                            Df Sum of Sq    RSS     AIC
## - density                   1     0.057 655.66 -1396.3
## - log_residual.sugar        1     0.150 655.76 -1396.0
## <none>                                  655.61 -1394.4
## - log_total.sulfur.dioxide  1     5.390 661.00 -1383.4
## - pH                        1     5.527 661.13 -1383.0
## - log_chlorides             1     7.075 662.68 -1379.3
## - volatile.acidity          1    34.832 690.44 -1314.0
## - log_sulphates             1    35.825 691.43 -1311.7
## - log_alcohol               1    72.703 728.31 -1229.0
## 
## Step:  AIC=-1396.26
## quality ~ volatile.acidity + log_residual.sugar + log_chlorides + 
##     log_total.sulfur.dioxide + pH + log_sulphates + log_alcohol
## 
##                            Df Sum of Sq    RSS     AIC
## - log_residual.sugar        1     0.094 655.76 -1398.0
## <none>                                  655.66 -1396.3
## - log_total.sulfur.dioxide  1     5.335 661.00 -1385.4
## - pH                        1     5.536 661.20 -1384.9
## - log_chlorides             1     7.147 662.81 -1381.0
## - volatile.acidity          1    34.812 690.48 -1315.9
## - log_sulphates             1    37.286 692.95 -1310.2
## - log_alcohol               1   108.388 764.05 -1154.7
## 
## Step:  AIC=-1398.03
## quality ~ volatile.acidity + log_chlorides + log_total.sulfur.dioxide + 
##     pH + log_sulphates + log_alcohol
## 
##                            Df Sum of Sq    RSS     AIC
## <none>                                  655.76 -1398.0
## - log_total.sulfur.dioxide  1     5.250 661.01 -1387.3
## - pH                        1     5.708 661.47 -1386.2
## - log_chlorides             1     7.056 662.81 -1383.0
## - volatile.acidity          1    34.720 690.48 -1317.9
## - log_sulphates             1    37.195 692.95 -1312.2
## - log_alcohol               1   112.828 768.59 -1147.3

## 
## Call:
## lm(formula = quality ~ volatile.acidity + log_chlorides + log_total.sulfur.dioxide + 
##     pH + log_sulphates + log_alcohol, data = dflo[c(2:10)])
## 
## Coefficients:
##              (Intercept)          volatile.acidity             log_chlorides  
##                  0.36690                  -0.93548                  -0.23784  
## log_total.sulfur.dioxide                        pH             log_sulphates  
##                 -0.08564                  -0.43301                   0.76832  
##              log_alcohol  
##                  3.10032

Информационный критерий AIC убирает 2 предикатора– density и log_residual.sugar.

Построим модель с использование минимального количества предикаторов.

lm.df.oecd <- lm.beta(lm(data = dflo[c(2:10)], quality ~ volatile.acidity + log_chlorides + log_total.sulfur.dioxide + pH + log_sulphates + log_alcohol))
summary(lm.df.oecd)

## 
## Call:
## lm(formula = quality ~ volatile.acidity + log_chlorides + log_total.sulfur.dioxide + 
##     pH + log_sulphates + log_alcohol, data = dflo[c(2:10)])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.61118 -0.35630 -0.04643  0.45756  1.90911 
## 
## Coefficients:
##                          Estimate Standardized Std. Error t value Pr(>|t|)    
## (Intercept)               0.36690      0.00000    0.50778   0.723 0.470065    
## volatile.acidity         -0.93548     -0.20592    0.10212  -9.161  < 2e-16 ***
## log_chlorides            -0.23784     -0.09472    0.05759  -4.130 3.82e-05 ***
## log_total.sulfur.dioxide -0.08564     -0.07463    0.02404  -3.562 0.000379 ***
## pH                       -0.43301     -0.08301    0.11658  -3.714 0.000211 ***
## log_sulphates             0.76832      0.21303    0.08103   9.482  < 2e-16 ***
## log_alcohol               3.10032      0.37947    0.18774  16.514  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6432 on 1585 degrees of freedom
## Multiple R-squared:  0.3616, Adjusted R-squared:  0.3592 
## F-statistic: 149.6 on 6 and 1585 DF,  p-value: < 2.2e-16

Все коэффициенты регрессии значимы (***), согласно модели quality будет больше при меньшем значении первых 4 признаков (volatile.acidity, log_chlorides, log_total.sulfur.dioxide, pH; максимальный по модулю– volatile.acidity) и большем значении 2 оставшихся (log_sulphates, log_alcohol; максимальный по модулю– log_alcohol).

Анализ остатков:

library(MASS)
library(car)
par(mfrow = c(2, 2))
plot(lm.df.oecd, 2) #Normal Q-Q
plot(lm.df.oecd, 4) #Cook`s distance
r <- stdres(lm.df.oecd)
jackknife.res <- studres(lm.df.oecd) #Residuals vs Deleted
plot(jackknife.res ~ r, main = "Residuals vs Deleted") + abline(coef = c(0,1))

## integer(0)

Outliers не выявлено (значения Cook`s distance не превосходит 0.05, на графике “Residuals vs Deleted” все точки лежат на прямой y=x), остатки распределены нормально, значит построенная модель хорошо описывает зависимость.

Спрогнозируем значение quality для 2 индивидов с использованием построенной модели линейной регрессии. Первый индивид будет обладать значениями признаков близкими к средним, второй соответствует близким к “идеальным” значениям признаков (для данной модели).

new_wine <- data.frame(volatile.acidity = 0.52, log_chlorides = log(0.08), log_total.sulfur.dioxide = log(46), pH = 3.311, log_sulphates = log(0.65), log_alcohol = log(10.42))
new_wine2 <- data.frame(volatile.acidity = 0.015, log_chlorides = log(0.015), log_total.sulfur.dioxide = log(10), pH = 3, log_sulphates = log(2), log_alcohol = log(14))
  
predict.lm(lm.df.oecd, new_wine, interval = "confidence")

##        fit      lwr      upr
## 1 5.654925 5.621168 5.688682

predict.lm(lm.df.oecd, new_wine, interval = "prediction")

##        fit      lwr      upr
## 1 5.654925 4.392829 6.917021

predict.lm(lm.df.oecd, new_wine2, interval = "confidence")

##    fit     lwr      upr
## 1 8.57 8.27884 8.861159

predict.lm(lm.df.oecd, new_wine2, interval = "prediction")

##    fit      lwr      upr
## 1 8.57 7.275195 9.864805

Как и ожидалось первый индивид получил оценку quality близкую к средней по выборке. А второй получил оценку и вовсе больше максимального по выборке (уж сильно “идеальными” характеристиками он обладает).